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We propose to combine the Trotter decomposition and a series expansion of the partition function 
for Hund's exchange coupling in a quantum Monte Carlo (QMC) algorithm for multiorbital systems 
that preserves spin and orbital rotational symmetries. This enables us to treat the Hund's (spin- 
flip and pair-hopping) terms, which is difficult in the conventional QMC method. To demonstrate 
this, we first apply the algorithm to study ferromagnetism in the two-orbital Hubbard model within 
the dynamical mean-field theory (DMFT). The result reveals that the preservation of the SU(2) 
symmetry in Hund's exchange is important, where the Curie temperature is grossly overestimated 
when the symmetry is degraded, as is often done, to Ising (Z2). We then calculate the tig spectral 
functions of S^RuCU by a three-band DMFT calculation with tight-binding parameters taken from 
the local density approximation with proper rotational symmetry. 

PACS numbers: 71.10.Fd; 71.15.-m; 71.27. +a 



I. INTRODUCTION 

As was pointed out by Slater as early as in the 1930s, 1 
the intra-atomic (Hund's) exchange interaction is of great 
importance for ferromagnetism in transition metals. Also 
the intensive study of transition-metal oxides — fueled by 
the discovery of high-temperature superconductivity — 
revealed the pivotal role of orbitals and Hund's ex- 
change, not only for ferromagnetism but, e.g., also for 
the metal-insulator transition 2-10 and unconventional 
superconductivity. 1 1-13 

For studying such correlation effects in multiorbital 
systems, the dynamical mean-field theory 14,15 (DMFT) 
is one of the most successful methods. It takes into ac- 
count temporal fluctuations but neglects spatial ones, 
and is known to describe the Mott-Hubbard transition. 15 
Recently, the multiorbital Mott-Hubbard transition has 
been studied most intensively by DMFT, 3-10 showing 
that the correct symmetry of the exchange coupling is es- 
sential, in particular, close to and at the Mott-Hubbard 
transition. 

Also for realistic material calculations, multiorbital 
DMFT is important. In particular, the so-called local 
density approximation (LDA)-fDMFT method 16 allows 
for calculating strong correlations in materials from first 
principles, and it has become one of the most active fields 
in solid state physics. 17 In this method, DMFT solves 
a model constructed from band structure obtained by 
the local density approximation. Since materials such 
as transition-metal oxides usually have several orbitals 
around Ep, the model is, in general, a multiorbital one. 
This multiorbital Hubbard model hence consists of the 
intra- and interorbital Coulomb interactions U and V , 
Hund's-exchange and the pair-hopping interactions J, as 
well as the one-electron hoppings obtained from the LDA. 

DMFT maps the model self-consistently onto an ef- 
fective impurity model which includes the same on-site 
interactions on an impurity site and an infinite number of 



bath sites, which are coupled to the impurity site through 
a hybridization (hopping). This impurity problem has to 
be solved numerically: For multiorbital models, the ex- 
act diagonalization 18 (ED) and the Hirsch-Fye quantum 
Monte Carlo (HFQMC) method 15 ' 19 are the most widely 
employed nonperturbative impurity solvers. Both ED 
and HFQMC methods have their pros and cons and are 
hence favorable in different situations. 

ED has the advantage that it can easily handle all the 
multiorbital interactions including Hund's coupling and 
the pair- hopping term. However, it is severely restricted 
in the number of bath sites, even in the Lanczos imple- 
mentations which can take into account most bath sites 
for T = or very low temperatures. 20 Since it cannot 
treat so many sites (and orbitals), DMFT (ED) studies 
have been usually limited to two-orbital systems. 

On the other hand, HFQMC simulation can handle not 
only two but more orbitals and, in contrast to ED, also 
produces continuous spectra. The former is important 
since there are various intriguing systems having three 
or more orbitals, especially ti g electron, perovskite-type 
transition-metal oxides, exemplified by the spin-triplet 
superconductor. S^RuOj 21 The continuous spectra are 
crucial for comparing with experiments, e.g., photoemis- 
sion spectroscopy. For these reasons the HFQMC method 
has been by far the most widely employed impurity solver 
as far as LDA+DMFT is concerned. 

In most DMFT HFQMC studies, 23 however, only the 
Ising (z) component of Hund's exchange coupling has 
been considered, with the pair-hopping term neglected. 
The reason is that the standard Hubbard-Stratonovich 
transformation 22 is only applicable to the density-density 
part (Hjj) of the interactions and something else is re- 
quired for the x and y components of Hund's coupling 
and the pair- hopping interaction {Hj). A straight- 
forward decoupling of Hj leads to a very severe sign 
problem, 2 making such QMC calculations virtually im- 
possible. 
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In order to take into account Hj 1 we previously pro- 
posed a transformation with a real and discrete auxiliary 
field for these terms in two-orbital systems. 12 ' 13 So far, s- 
wave superconductivity 13 and the orbital-selective Mott 
transition 9 ' 10 have been investigated with our transfor- 
mation. 

Despite these successes, some problems remain to be 
solved. First, it is still difficult to treat more than two 
orbitals: Since the exponential of Hj is not equal to the 
product of the exponentials of the distinct two-orbital 
parts Hj lm (m and m' denote orbitals), 25 we cannot 
decouple Hj for three or more orbitals by the simple 
product of two-orbital parts. Second, the fermionic sign 
problem 24 caused by Hj prevents us from studying low 
temperatures. 

In the present paper we propose and show that these 
difficulties are greatly remedied by employing a pertur- 
bation series expansion (PSE), separating out the prob- 
lematic Hj term. It should be noted that the method 
is numerically exact (and nonperturbative) . Also, the 
method is easily implemented, since our updating algo- 
rithm for sampling weights and Green's function is ba- 
sically the same as the Hirsch-Fye algorithm, with the 
differences being only for a factor multipling the weights 
and an additional value for the auxiliary field. While we 
focus here on DMFT, i.e., QMC simulations for impurity 
problems, the same ideas can be employed also for the 
lattice QMC method. 

The auxiliary-field QMC method based on PSE for 
electron systems was first proposed by Rombouts et al. 26 
These authors applied PSE to the finite-size single-orbital 
Hubbard model and succeeded in obtaining results 
without time-discretization errors with less computa- 
tional time than the conventional, Trotter-decomposition 
algorithm 7 bs81 Although the method uses PSE, it 
counts all the contributions of the interaction since the 
perturbation order taken into account is higher than the 
maximum order of the samples above which the weight 
is virtually zero. So the scheme is essentially nonper- 
turbative. Rubtsov et al. 2S proposed another algorithm 
to evaluate the series expansion of the partition function 
and applied it to solve the DMFT equations. It does 
not involve any auxiliary fields but uses Wick's theorem. 
Recently Werner et al. 29 proposed to use a perturbation 
expansion starting from the strong-coupling limit. 

In previous work, 30 we extended Rombouts' algorithm 
to the multiorbital Hubbard Hamiltonian, using a similar 
transformation as in Rcf. 13 [Eq. (3)] for Hj, i.e., we 
expanded the Boltzmann factor (operator) with respect 
to the total interaction Hjj + Hj shifted by a constant. 
Then, we discretized the imaginary time (3 = LAt and 
used a Hirsch-Fye- like updating algorithm for solving the 
impurity model in the DMFT context. Although the 
method significantly relaxes the sign problem and allows, 
in principle, for more than two orbitals, it turned out that 
the calculations are too heavy at low temperatures or for 
strong couplings, especially for more than two orbitals. 
That is because the computational effort increases with 



the perturbation orders appearing in the Monte Carlo 
samples. This order can become very large (see Fig. 2) 
in multiorbital systems since there are many interactions: 
(2M - 1)M terms in H v and 2(M - 1)M terms in Hj 
per site, where M is the number of orbitals. 

To overcome this difficulty, we here propose to combine 
the HFQMC and the PSEQMC methods, i.e., to adopt 
the series expansion for Hj, while the standard Trotter 
decomposition and Hubbard-Stratonovich decoupling are 
employed for Hij. This algorithm enables us not only to 
handle three or more orbitals but also to reach much 
lower temperatures or stronger couplings than HFQMC 
(Rcf. 13) or PSEQMC calculations, 30 even for the two- 
orbital Hubbard model. 

In the following, we first introduce the multiorbital 
Hubbard model in Sec. II to turn to our algorithm in 
Sec. Ill, enclosing details for the weight factor in the 
Appendix. A validation and benchmarks are presented 
in Sec. IV. Applications to local spin moments and to 
ferromagnctism in the two-orbital Hubbard model along 
with an application to Sr 2 RuC>4 are discussed in Sec. V. 
Section VI summarizes the paper and gives an outlook. 



II. MODEL 

The M-orbital Hubbard Hamiltonian reads 

H =H a + Hu + Hj, (1) 
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where cj m(T (Q m(T ) creates (annihilates) an electron with 
spin d in orbital tyi at site i, and nima = c\< ma Cimo- 
describes the hopping of electrons between two neigh- 
boring sites (ij), which we assume, unless otherwise in- 
dicated, to be orbital independent. Hu is defined to be 
all the density-density interactions, i.e., the intraorbital 
(U) and interorbital ([/') Coulomb interactions and the z 
(or Ising) component of Hund's coupling J. By contrast, 
Hj takes care of the terms that cannot be written as 
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density-density interactions, which consist of the x and y 
component of Hund's coupling as well as the pair- hopping 
interaction (second term), in which two electrons in an 
orbital transfer to other orbitals. The Hamiltonian is 
rotationally invariant not only in spin space but also in 
real (orbital) space if we satisfy the condition U = U'+2J, 
which comes from the symmetry preserved between the 
Coulombic matrix elements for orbitals in a central field. 
We postulate the rotational invariance henceforth. 

In the DMFT approximation, the Hamiltonian (1) is 
approximated by an impurity problem where the inter- 
actions Hu and Hj are restricted to a test atom (the 
impurity) embedded in the systems; and Hq describes 
the coupling to a noninteracting bath which has to be 
determined self-consistently. 31 When we develop our al- 
gorithm in the next section, Hq, Hj, and Hjj denote 
these terms of the impurity model, but the same ideas 
can also be applied for lattice QMC simulations of the 
Hubbard model. 



III. (TROTTER + SERIES-EXPANSION) 
METHOD 



We start with the series expansion of the Boltzmann 
factor for continuous time after Ref. 26. However, here 
we perform this only for Hj, i.e., 

e f-f3H = e -/3(H +HL0+7-/3H, 
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where we have shifted the Boltzmann factor by a constant 
7 for (3Hj to apply the auxiliary-field transformation Eq. 
(9) below. 

Now we discretize the imaginary-time integrals and 
with the notation Xi = 7 — (3Hj, Eq. (2) equals 



if L is taken to be greater than the maximum perturba- 
tion order k max (defined and displayed below) appearing 
in the Monte Carlo samples, so that there are no con- 
tributions from higher-order terms. In practice, we can 
make L much larger than k max (see Fig. 2 below): k max 
depends on Hund's coupling J, where J is physically not 
so large, whereas we can choose L to satisfy L > (3U. 

Second, we replace those terms having consecu- 
tive Xi's in Eq. (3) by proximate terms includ- 
ing only one Xi's per imaginary time interval At. 
For example, • • • XiXie~ AT ( Ho+Hu ) ■ ■ ■ is replaced by 
• • • X 1 e^ AT< -"<' + " u ^X 1 ■■■ . This replacement reduces the 
number of possible configurations remarkably and casts 
the summation (3) into the form (4) similar to the Trot- 
ter decomposition, which enables us to employ their stan- 
dard Hirsch-Fye algorithm with only a slightly more com- 
plicated auxiliary field at each time slice. The error of 
this approximation (commutation) is O(Ar), i.e., of the 
same order as the time discretization, as long as the aver- 
age order of the scries expansion (k) is sufficiently smaller 
than L. This is simply because the terms having two or 
more consecutive Xi's rarely appear for (k) <C L. For 
example, consider the second order terms in Eq. (3). Al- 
together, there are L(L + l)/2 second-order terms, but 
only L of these terms have two consecutive Xi's with 
the same imaginary time interval. Hence, the error is 
at most 0(2At/L). Similar argument for higher orders 
justify the replacement as long as (k) <C L. Since we 
do not simply drop the terms with two or more consec- 
utive Xi's, but replace them by terms where the Xi's 
are shifted to neighboring imaginary time intervals, we 
have to multiply the Boltzmann factor by a factor F to 
account for these replacements. The detailed derivation 
of F is given in the Appendix. 

Now, we separate out Hjj in Eq. (4) using the Trotter 
decomposition as 



e -Ar(H +Hu) = e -ArH„ e -ArHu + {At 2 ), (5) 
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We now show that this sum can be rewritten as 
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where F is a positive weight factor, k = s i> an d 

Xq = 1. To obtain the representation (4), we first cut off 
the k summation in Eq. (3) at L. This cutoff is justified 



so that Eq. (4) has a similar form to the standard 
HFQMC method. The e - ArHu term is then decoupled, 
as usual, into a sum of one-body exponentials with the 
Hubbard-Stratonovich transformation, 
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where V stands for U, U', or U' - J, X v = ln( e l Ary l/ 2 + 
\/el Ary l — 1). We have also displayed the case of at- 
tractive interaction, which we shall require when we do 
the procedure described in Ref. 33. Including all the 
(2M — l)M interactions of density-density type, the de- 
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coupling for e ^ tHu is given by 
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where P(= 1,...,N V = 2^ 2M ^ M ) designates configu- 
rations of the auxiliary-field set ({p m }, {q™ m }, {r™ m }) 
with p m , q™ m , r™ m = ±1 denoting the fields for the U, 
V, and U 1 — J terms, respectively. 

For X\ = "f—f3Hj in Eq. (4) we construct an auxiliary- 
field transformation as follows. 30 ' 32 We first decompose 
Xi into the sum of all distinct two-orbital parts as 

7 - pAj = [7™" 

m<m' 

We then apply the decoupling 
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to every pair of orbitals, where 



c \na C m'a + C^Cnw, 
1, 1 + K 

- In , 

2 f — K 



(3J 



ry'l 



J < 1. 



Note that the difficulty in the HFQMC calculations com- 
ing from the noncommutativity of H 7 p m 's is lifted in Eq. 
(8), so that we can readily deal with more than two or- 
bitals. Combining Eqs. (8) and (9), we end up with 
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where 5[= f , A/ = 4M(M— f )] corresponds to the set 
(s,tf,ti) for all the M(M— 1)/2 pairs of orbitals (m, m'). 

Collecting the addends from the decoupled Hjj and Hj 
terms, we finally obtain 
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FIG. 1: (Color online) Real and imaginary parts of the self- 
energy against the Matsubara frequency w„ for the two-orbital 
Hubbard model for (a) an insulating case with n = 2, f3 = 
10, U' = 2, J = 0.4, and (b) a metallic case with n = l,/3 = 
6, 17' = 2, J = 1. Result with the Hirsch-Fye algorithm (Ref. 
13) is shown in black squares, and the present QMC result in 
red triangles. 



with Qq = 1 and Sj = (for Si = 0) or 1 (otherwise). 33 
Note that because F(0;0, ...,0) = 1, the zeroth-order 
term in Eq. (11) reproduces the Hirsch-Fye algorithm 
with Ising-type Hund's coupling. 

Owing to the form of Eq. (11), we can apply the same 
algorithm as in HFQMC for Monte Carlo sampling. Even 
the updating equations for single spin flips are the same. 



IV. BENCHMARKS 

As a benchmark, we compare in Fig. 1 the electron 
self-energy obtained from our algorithm with that from 
the HFQMC method of Ref. 13 for the two-orbital Hub- 
bard model. In these QMC methods the self-energy is 
obtained from Dyson's equation, where Green's function 
is calculated by averaging over QMC samples. We chose 
the hypercubic lattice with bandwidth W — 2, and took 
6 x 10 6 Monte Carlo samples for both methods. We can 
see that the two results agree with each other within error 
bars for both (a) an insulating case at half filling n = 2 
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FIG. 2: (Color online) The distribution N k of the order 
of perturbation k for the two-orbital Hubbard model with 
n = 1.9, j3 = 8,U' = A, J = 0.2 obtained with the present 
algorithm (L = 64 or 100) and with the PSEQMC algorithm 
(L = 150) (Ref. 30). 



FIG. 3: (Color online) Computable regions (hatched) in the 
T-J parameter space that can be computed with the present 
and with the Hirsch-Fye QMC method (Ref. 13) for the two- 
orbital Hubbard model with U' = 4, W = 2. Here we define 
the computability by requiring average signs to be greater 
than 0.01. 



with f3 = 10, U' = 2, J = 0.4, L = 100, and (b) a metallic 
case at n = l,p = 6,V = 2, J = 1,L = 64. We no- 
tice, however, that the statistical error is much smaller 
in the present QMC approach than in the HFQMC one. 
This is because the number of negative signs is greatly 
reduced in the present scheme: The sign problem is mit- 
igated. Quantitatively, the average sign in the QMC 
weights is 0.01 (0.03) for HFQMC methods while they 
are increased to 0.30 (0.50) in the present algorithm in 
case (a) [(b)]. This also implies that the present method 
can reach much lower temperatures. We note that, while 
7 is arbitrary, we adopted in these and following calcula- 
tions 7 — (3V ~ 0.1 — 0.3, which has turned out to reduce 
the computational time to some extent. 

Figure 2 depicts a typical distribution Nk of the or- 
der of perturbation k contributing in the Monte Carlo 
simulation for n = 1.9, f3 = 8, U' = 4, J = 0.2. We 
can immediately see that the present algorithm has a 
peak in the distribution residing at much lower k than 
for the PSEQMC approach. 30 This is natural, since the 
present method uses the expansion only with respect to 
Hj, while the PSEQMC method expands with respect to 
the total interaction Hu + Hj. We notice that the weight 
is virtually zero above a certain order, k max , in actual 
calculations, so that, although the method is called per- 
turbation expansion, it takes account of all orders in fact. 
We find that the maximum order in the distribution is 
k max ~ 100 for the PSEQMC method, while k ma x ~ 15 is 
much lowered for the present QMC method. This means 
that L must be taken to be > 100 for the PSEQMC 
method to take care of all orders, while L > (3U ~ 35 suf- 
fices for the present algorithm to take care of all orders. 
Such a smaller value of L dramatically reduces the com- 
putational effort in QMC simulations, which increases 
proportionately to L 3 . Moreover, the average order (k) 
is about 4 for the present QMC, which means that the 
approximation employed to obtain the form (4) has only 
a very minor effect on the results, and hence is justifiable. 



This can also be confirmed from the fact that the results 
do not significantly depend on L. 

Figure 3 shows the computable regions for the present 
QMC and for the HFQMC methods (Ref. 13) when Hj is 
included. Here we define the region as computable when 
the average sign is greater than 0.01. We can see that 
a much wider parameter region becomes computable in 
the present algorithm than in the HFQMC method. For 
small J (< 0.2), we can explore five to ten times lower 
temperatures. We attribute this improvement to the fact 
that Hj (which is the source of negative weights) appears 
L times for every sample in the HFQMC method, while 
we have only (k) such terms on average in the present 
QMC algorithm. 



V. APPLICATIONS 



A. Local spin moment 

Let us now turn to the first application of our algo- 
rithm. We first examine how results for the Hamiltonian 
with and without Hj, i.e., with Z2 and SU(2) symmetry, 
respectively, can differ. In Fig. 4, we show the tempera- 
ture dependence of the local spin moment in the param- 
agnetic phase at half filling in the two-orbital Hubbard 
model for U' = 3, J = 0.2 on a hypercubic lattice with 
bandwidth W = 2. At these parameters, the system 
becomes insulating at zero temperature for both the Z2 
(Ising) and the SU(2) cases. The spin moments are de- 
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FIG. 4: (Color online) Local spin moment for the two-orbital 
Hubbard model with SU(2)-type Hund's and pair-hopping in- 
teractions as compared with Ising-type Hund's coupling for 
U' = 3, J = 0.2, n = 2 (half filling) and the Gaussian den- 
sity of states of width W — 2. Solid lines are guides to the 
eye, dashed lines extrapolations, and the dotted line the local 
moment for free electrons. 
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For the Ising case (S%) is always larger than (S%. y ), 
while for the SU(2) case these components have the same 
value (within statistical error bars) as they should, which 
lies between the (S z ) and (S x ), (S y ) in the Ising case. As 
the temperature is lowered in the Ising case, (S z ) and 
(S%. y ) approach 1 and |, respectively, while for SU(2) 
symmetry both expectation values converge to | . This is 
easily understood by considering the fact that the ground 
state in the atomic limit is a doublet with S z — ±1 in 
the Z2 case, whereas it is a triplet with S z — 0, ±1 in the 
SU(2) case. 

B. Ferromagnetism 

We now turn to the long-range ferromagnetic order- 
ing, which reveals a more important difference between 
Z 2 and SU(2) symmetries as Fig. 5 shows. To this end, we 
have calculated the temperature dependence of the mag- 
netic susceptibility in the two-orbital Hubbard model for 
U' = 2.5, J = 1, n = 1.25 with a semielliptical density of 
states of width W — 2. The lattice susceptibilities are ob- 
tained here through the Bethe-Salpeter equation. 15 The 
Curie temperature for the Hamiltonian with Ising-type 
Hund's coupling is estimated by a Curie- Weiss extrapo- 
lation to be 0.01-0.02, in accordance with the previous 
results. 24 In contrast, the present result shows that the 
inverse susceptibility does not intersect the horizontal 



FIG. 5: (Color online) Inverse spin susceptibility for the two- 
orbital Hubbard model with SU(2)-type Hund's and pair- 
hopping interactions as compared with the Ising-type Hund's 
coupling for U' = 2.5, J = 1 with the semielliptical density of 
states of width W — 2. The solid lines are guides to the eye, 
and the dashed lines extrapolations. 
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FIG. 6: (Color online) Quasiparticle density of states for ti a 
orbitals in Sr2Ru04, obtained by a simplified LDA+DMFT 
(QMC) calculation with three-orbital rotationally invariant 
Hamiltonian (1). 



axis, implying an absence of ferromagnetic transitions, 
for these parameters if the SU(2) symmetry of Hund's 
exchange and the pair-hopping term are respected. 

The result suggests that an Ising-type treatment of 
Hund's coupling grossly overestimates the tendencies to- 
ward ferromagnetic ordering. Indeed, Ising-type DMFT 
calculations for manganites 34 and iron 35 gave Curie tem- 
peratures much higher than the experimental results. 
The present comparison indicates that the main reason 
for overestimating the Curie temperature in such calcu- 
lations is not only the mean-field nature of DMFT but 
the incorrect symmetry of the Hund's exchange. For an 
Ising-type Hund's exchange, (local) transverse spin fluc- 
tuations, which are the main source of softening magnetic 
moments, are suppressed. 



C. Quasiparticle density of states for Sr2Ru04 

Finally we demonstrate that the present QMC method 
is actually applicable to realistic three-orbital systems. 
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As an illustration, we take Sr2Ru04, which can be 
treated as a (t^g) three-orbital system, consisting of a 
wide two-dimensional band (d xy ) and two degenerate 
narrow one-dimensional bands (d xz ,d yz ). The material 
belongs to an important family of transition-metal oxides 
where superconductivity, 21 magnetism, 37 and orbital- 
selective Mott transition 4-8 are discussed. To our knowl- 
edge, this is the first QMC calculation for a three-orbital 
system with SU(2)-symmetric Hund's exchange and pair- 
hopping interaction. 

Licbsch and Lichtenstein 36 already obtained the quasi- 
particle density of states for this material with a simpli- 
fied LDA+DMFT (HFQMC) calculation using the three- 
orbital Hamiltonian without Hj (i.e., an Ising-type treat- 
ment of Hund's coupling). As a LDA input they con- 
structed a tight-binding model that reproduces the LDA 
density of states. We have employed the same initial 
tight-binding model (for the density of states), the same 
number of t2 g electrons per site (n = 3.7), and the same 
interaction parameters V = 0.8 cV, J = 0.2 cV, but a 
somewhat higher temperature T — 25 meV. 

The LDA+DMFT spectral functions, obtained with 
the maximum entropy method, are shown in Fig. 6. The 
quasi-two-dimensional (2D) d xy band has a van Hove sin- 
gularity just above the Fermi energy Ep, while the quasi- 
1D d xz and d yz bands have broad peaks on both sides of 
Ef and a small peak near Ef- Although these structures 
are similar to those of Liebsch and Lichtenstein, 36 it may 
be due to the small value of J employed here. In gen- 
eral, Hj can have a significant influence on quasiparticle 
spectra. We have obtained results substantially changed 
by Hj, for example, enhancement of Kondo resonance 
by Hj, for artificially elevated U, U 1 , and J (not shown 
here). Actually, values of U , U' , and J larger than those 
used here, as usually assumed for transition-metal ox- 
ides, were recently suggested for Sr 2 RuC>4, 38 and it will 
enhance the effect of Hj. 



VI. SUMMARY AND OUTLOOK 

In summary, we have formulated a numerically exact 
auxiliary-field QMC scheme for the spin- and orbital- 
rotationally-invariant Hamiltonian (1) by combining se- 
ries expansion and Trotter decomposition. We have im- 
plemented this algorithm for impurity models and em- 
ployed it to solve the DMFT equations. The approach 
enables us not only to address three- (or more-)orbital 
systems but also to reach much lower temperatures than 
the previous QMC methods, which is the case with two- 
orbital systems as well. 

The calculation of the magnetic susceptibility shows 
that an Ising treatment of Hund's coupling overestimates 
tendencies towards ferromagnctism, which we attribute 
to neglected transversal spin fluctuations. In this sense 
the SU(2) symmetry in Hund's exchange is indicated to 
be mandatory for an accurate study of magnetic phenom- 
ena with (LDA+)DMFT. We finally applied the present 



method to a realistic three-orbital system Sr2Ru04, to 
demonstrate that the QMC algorithm can be employed 
for LDA+DMFT calculations. 

Since we can apply the usual Hirsch-Fye updating algo- 
rithm with only minor changes in the additional values 
of the auxiliary field and a factor for the weights, the 
present scheme can be conveniently applied to multior- 
bital DMFT studies, especially for LDA+DMFT calcu- 
lations requiring three orbitals. Although the fermionic 
sign problem still remains at low temperatures, we can 
reach much lower temperatures than before. 13 Also, a 
combination of the present method with the projective 
QMC method 10 ' 39 is very promising, which may provide 
an avenue for examining ground states, which is now un- 
der way. 
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Appendix 

We derive here the factor F(k; s\, S2, Sl) (si — 
0,1; k = 5Z j=1 s i) f° r Eq. (4). We introduce this fac- 
tor to account for the contribution from the terms with 
consecutive AVs at the same imaginary-time interval in 
the sum (3). These terms have been replaced with terms 
having X\ 's on proximate imaginary-time intervals in Eq. 
(4). In the following we abbreviate e ~ AT ( H o+Hu) as ^ 
and Xi as x. 

Central to our considerations are those terms that 
include a substring xhxhx ■ ■ ■ hx where x and h ap- 
pear alternately m and m — 1 times, respectively, with 
2 < m < L. This is because any term with consecutive 
x's will be replaced with a term having such substrings. 
For example, xxhhx ■ ■ ■ hx and xhxxh ■ ■ ■ hx are both ap- 
proximately the same as xhxhx. ..hx [where commuting 
an x and a h yields an error <~ O(Ar)]. In general, we 
commute x's and h's until there are no consecutive x's 
any longer. Hence, we end up with an alternation of x's 
and /i's, i.e., an xhxhx... hx substring. Because of these 
replacements, terms having such a substring have to be 
weighted more. In the following, we construct a rule for 
the replacement and weighting factor, avoiding a double 
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counting. 

Let us denote by i a position (from the left) in a sub- 
string which consists of m x's and (m— 1) /i's altogether. 
We define Qj and as the number of x's and /i's that 
are in [l,i]. All substrings having 

on > Pi for all i (15) 

will be replaced with xhxhx...hx which has alternating 
x's and /i's. 

For example, for m = 2 xhx and xxh will be replaced 
with xhx; for m = 3, xhxhx,xhxxh,xxhhx,xxhxh and 
xxxhh will all be replaced with xhxhx. The condition 
(15) is necessary to avoid a double counting. However, 
the condition (15) excludes the substrings situated at the 
end of the imaginary time interval for which 3 i', so that 
< Pi'. These terms are replaced with a substring 
xhxhx ■ ■ ■ hx, where the last x is at L, namely, only when 
the last x is at L, does xhxhx ■ ■ ■ hx replace all the sub- 
strings having m x's and (m — 1) /i's, which requires a 
separate treatment (c factors below). 

A second factor to be taken into account is the cor- 
rection of the volume in the imaginary-time integrals; 
namely, the weight for those terms having consecutive 
x's at 

k 7^ = Ji+2 = ■■■ = ji+i ^ ji+i+i (16) 

in the sum (3) should be reduced by a factor jf, since the 
imaginary times originally satisfy a relation 

U+i < t i+ 2 < ••• < U+i (17) 

in Eq. (2). Hence the volume for I consecutive x terms, 
as in Eq. (16), should be reduced to 

Let us now introduce the quantity b(i,j) for the weight 
(apart from the volume factor L~i) of all the substrings 
having i ft's and j x's. For j consecutive x's, we simply 
have the aforementioned 1/1! factor, i.e., 

6(0,j) = i for 0<j<L. (18) 

This is the starting point for the recurrence formula, 

j 1 

^'■?) = E 7^— Mi 6 (i-1.*) for 1 < * < j < L, (19) 



which arises from taking away the rightmost elements 
of the type hxxx ■ ■ ■ x with (j — k) x's from the sub- 
string of length i+j. At the end, the recurrance formula 
Eq. (19) yields the weighting factor at for the substring 
l xhxhx ■ ■ ■ hx' with m x's and m — 1 h's: 

di = b(i-l,i) = b{i,i) for 1 < i < L. (20) 

Only when the last x in xhxhx... hx is situated at the 
end of the imaginary-time interval (L) do we use the fac- 
tor c m instead, which is obtained via the slightly different 
recurrence formula [see Fig. 7(b)], 

d(0,j) = j for 0<j<L, 

j 1 

k=o u >' 
for 1 < i < L - 1 and < j < L, 

Ci = d(i for 1 < i < L. (21) 

From the a's and c's, the total weight F is calculated 
by multiplying the contributions a m and c m from each 
xhxhx... hx-type substring in the Boltzmann factor and 
the volume L~ k . For example, for L = 8, 

J F 1 (2;l,l,0,0,0,0,0,0) = a 2 L- 2 , 
F(5;0,l,0,l,0,l,l,l) = csL- B , 
F(6;l,0,l,l,0,l,l,l) = O2C 3 L- 6 . (22) 

In the first example, the Boltzmann factor is 
hxhxhhhhhh. This array replaces itself and 
hxxhhhhhhh, which is weighted ^j, and thereby 
the factor is a 2 = 1 + ^ = § multiplied by L~ 2 . The 
second example corresponds to hhxhhxhhxhxhx, where 
the substring to be multiplied by a factor is only the 
last part xhxhx. Therefore the factor is C3-L~ 5 . In the 
last example hxhhxhxhhxhxhx, two substrings xhx and 
xhxhx are multiplied by a 2 and C3, respectively. So the 
total factor is a 2 x c^L^ 6 . 
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